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Abstract. Structure-preserving numerical schemes for a nonlinear parabolic fourth- 
order equation, modeling the electron transport in quantum semiconductors, with peri- 
odic boundary conditions are analyzed. First, a two-step backward differentiation formula 
(BDF) semi-discretization in time is investigated. The scheme preserves the nonncgativity 
of the solution, is entropy stable and dissipates a modified entropy functional. The exis- 
tence of a weak semi-discrete solution and, in a particular case, its temporal second-order 
convergence to the continuous solution is proved. The proofs employ an algebraic relation 
which implies the G-stability of the two-step BDF. Second, an implicit Euler and g-step 
BDF discrete variational derivative method are considered. This scheme, which exploits 
the variational structure of the equation, dissipates the discrete Fisher information (or 
energy). Numerical experiments show that the discrete (relative) entropies and Fisher 
information decay even exponentially fast to zero. 



1. Introduction 

This paper is devoted to the study of novel structure-preserving temporal higher-order 
numerical schemes for the fourth-order quantum diffusion equation 



(1) n t + div J^nV (~^-J J = °> ^^t>0, n(0) = n Q , 

where T d is the d- dimensional torus. This equation is the zero-temperature and zero-field 
limit of the quantum drift-diffusion model, which describes the evolution of the electron 
density n(t) = n(t,-) in a quantum semiconductor device; see [21]. It was derived in [7] 
from a relaxation-time Wigner equation using a Chapman-Enskog expansion around the 
quantum equilibrium. For smooth positive solutions, (CD) can be written in a symmetric 
form for the variable logn: 

(2) n t + -<9j(n<9jlogn) = 0, x e T d , t > 0, n(0) = n , 



Date: August 28, 2012. 

2000 Mathematics Subject Classification. 65M06, 65M12, 65M15, 35Q40, 82D37. 

Key words and phrases. Derrida-Lebowitz-Spcer-Spohn equation, discrete entropy-dissipation inequa- 
lity, Fisher information, BDF time discretization, numerical convergence, discrete variational derivative 
method. 

The first and last author acknowledge partial support from the Austrian Science Fund (FWF), grants 
P20214, P22108, and 1395, and the Austrian-French Project of the Austrian Exchange Service (OAD). 

1 



2 



MARIO BUKAL, ETIENNE EMMRICH, AND ANSGAR JUNGEL 



where here and in the following, we employ the summation convention over repeated indices 
and the notation di = d/dxi, dfj = d 2 /dxidxj. This is the multidimensional form of the 
so-called Derrida-Lebowitz-Speer-Spohn (DLSS) equation. Its one-dimensional version was 
derived in [8] in a suitable scaling limit from the time-discrete Toom model and the variable 
n is related to a limit random variable. 

The main difficulties in the analysis of ([I]) (or fl2])) are the highly nonlinear structure, 
originating from the quantum potential term /A^/nj ^/n in (pQ), and the fourth-order differ- 
ential operator, which lacks a maximum principle. 

These difficulties have been overcome by exploiting the rich mathematical structure of 
([2]). First, equation (T5]) preserves the nonnegativity of the solutions [22]: Starting from 
a nonnegative initial datum, the weak solution stays nonnegative for all time. Second, 
(T5]) allows for a class of Lyapunov functionals and so-called entropy dissipation estimates. 
More precisely, the functionals 

E a [n} = —- -/ n a dx (a ^0,1), E x \n] = / (raflogra - 1) + l)dz 

a(a - 1) J T d Jjd 

are Lyapunov functionals along solutions to (j2J), i.e. dE a [n]/dt < if (y/d — l) 2 /(d + 2) < 
a < (y/d + l) 2 /(d + 2), and the entropy dissipation inequality 

(3) -r,E a [n] + 2n a [ (An a/2 ) 2 dx < 

dt Jjd 

holds if (y/d - l) 2 /(d + 2) < a < (y/d + l) 2 /(d + 2). The constant n a > can be 
computed explicity, see Lemma |6] below. For a — 1, inequality ()3]) can be interpreted as 
the dissipation of the physical entropy. Third, equation (JTJ) is the gradient flow of the 
Fisher information 

(4) F[n) = [ |Vv^| 2 dx 

with respect to the Wasserstein metric [H]. As the variational derivative of the Fisher 
information equals 5F[n]/5n = — Ay^t/^/n, a straightforward computation shows that the 
Fisher information is dissipated along solutions to ([TJ, 

(5) ±F[n] + [ n v(^) "dx = 0. 

dt J T d \ on / 

Since the Fisher information can be interpreted as the quantum energy, the latter can be 
seen as an energy dissipation identity. 

Whereas the local-in-time existence of positive classical solutions for strictly positive 
W 1,P (T ) initial data with p > d could be proved using semigroup theory [2], global-in- 
time existence results were based on estimates ((3]) and (jSJ). More precisely, the global 
existence of a nonnegative weak solution was achieved in [21] in the one-dimensional case. 
This result was extended later to several space dimensions in [22], employing entropy 
dissipation inequalites, and in [T4] , exploring the variational structure of the equation. 

From a numerical viewpoint, it is desirable to design numerical approximations which 
preserve the above structural properties like positivity preservation, entropy stability, and 
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entropy or energy dissipation on a discrete level. For a constant time step size r > 0, 
let tk = kr (k > 0). If approximates the solution n(tk) to §2§ at time t^, we call a 
numerical scheme entropy dissipating if E a [nk+i] < E a [rik] for all k > with a in a certain 
parameter range, and entropy stable if there exists a constant C > such that E a [rik] < C 
for all k > 0. In this paper, we investigate the entropy stability and entropy dissipation of 
backward differentiation formulas (BDF). 

In the literature, most of the numerical schemes proposed for ([2]) are based on an im- 
plicit Euler discretization in one space dimension. In [25] , the convergence of a positivity- 
preserving semi-discrete Euler scheme was shown. A fully discrete finite-difference scheme 
which preserves the positivity, mass, and physical entropy was derived in [I] . During et al. 
[H] employed the variational structure of ([2]) on a fully discrete level and introduced a dis- 
crete minimizing movement scheme. This approach implies the decay of the discrete Fisher 
information and the nonnegativity of the discrete solutions. Finally, a positivity-preserving 
finite-volume scheme in several space dimensions for a stationary quantum drift-diffusion 
model was suggested in [5]. 

Positivity preserving and entropy consistent numerical schemes have been investigated in 
the literature also for other nonlinear fourth- and second-order equations. For instance, a 
positivity preserving finite difference approximation of the thin-film equation was proposed 
by Zhornitskaya and Bertozzi [21]. Finite element techniques for the same equation were 
employed by Barrett, Blowley, and Garcke [T], imposing the nonnegativity property as a 
constraint such that at each time level a variational inequality has to be solved. Further- 
more, entropy consistent finite volume-finite element schemes were suggested and analyzed 
by Griin and Rumpf [T7J [18]. Furihata and Matsuo [13] developed the discrete variational 
derivative method to derive conservative or dissipative schemes for a variety of evolution 
equations possessing a variational structure. Entropy dissipative fully discrete schemes for 
electro- reaction-diffusion systems were derived by Glitzky and Gartner [16] . 

In most of these works, the time discretization is restricted to the implicit Euler method, 
motivated by the fact that the solutions often lack regularity. However, high-order schemes 
often still yield smaller time errors than the Euler scheme, and this improved accuracy is 
vital to match the spatial approximation errors. A difficulty of the analysis is that the time 
discretization has to be compatible with the entropy structure of the equation. This is the 
case for the first-order implicit Euler discretization. Indeed, multiplying the semi-discrete 
scheme 



where r > is the time step and approximates n{tk) with t^ = rk, by log rifc +1 and 
using the elementary inequality 



(which follows from the convexity of x H- xlogx), it was shown in j22j Lemma 4.1] that 





(7) 






k > 0. 
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As a consequence, k h- > E a [n k ] is nonincreasing and the entropy dissipation structure is 
preserved. It is less clear whether higher-order approximations yield entropy dissipating 
numerical schemes. In this paper, we prove this property for the two-step BDF method. 

Two-step BDF (or BDF2) methods have been employed in the literature to approximate 
various evolution equations in different contexts. We just mention numerical schemes for 
incompressible Navier-Stokes problems [101 [151 HSj , semilinear and quasilinear parabolic 
equations [ITJ [28] , and nonlinear evolution problems governed by monotone operators [T2l 
[20] . To our knowledge, temporal higher-order schemes for the quantum diffusion equation 
([T|) have been not considered so far. 

In the following, we detail our main results. First, we analyze the BDF2 time approxi- 
mation of the DLSS equation, written in the form 

(8) l n W2 K / 2)t + l d %{ndl logn) = 0, 

which was already used in [27] in a different context. Introducing the variable v k '■= n^ 2 , 
which approximates n(t k ) a ^ 2 , the semi-discrete BDF2 scheme for ([8]) reads as 

(9) J^T 1 U^+i - 2«* + + \dl(n k+1 dl \ogn k+1 ) = in T d , k > 1. 

Here, vq = n^ 2 is given by the initial datum n , and v\ is the solution to the implicit Euler 
scheme 

(10) — v 2 1 /a - 1 (v 1 -v ) + l-d* J (n 1 dl\ogn l ) = in T d . 

OCT 2 

The existence of a weak solution to the scheme (!9|)- (fT0|) is provided by the following 
theorem. 

Theorem 1 (Existence of solutions and entropy stability). Let l<<i<3 ; l<a< 
(Vd + l) 2 /(d + 2), and let no e L 3 (T d ) be a nonnegative function. Then there exists a 
weak solution v\ = n^ 2 of the implicit Euler scheme ( flOl) and a sequence (vk) = { n t^ 2 ) 
of weak nonnegative solutions to satisfying v k > in T d , v k G H 2 (Y d ), and for all 

<p e W 2 >°°{T d ), 

(11) — / vt+i" 1 - 2v k + 7; v k-i) <Pdx 

J T d \2 2 J 

If a > 1, the scheme is entropy stable and the a priori estimate 

4 m f 

(12) E a [n m ] + -K a rV / (A« /2 )) 2 dx < E Q [n ], m > 1, 
holds, where K a > is defined in Lemma\^ 
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When we redefine the entropy, we are able to prove entropy dissipation of the semi- 
discrete scheme. For this, introduce the modified entropy 

ES[n H ,n k - 1 ]= 1 / {nt + (2nf-n a k / _\) 2 )dx, k > 1. 
2a(a - 1J J Jd 

This definition is motivated by the inequality 

2 ( -a - 2b + -c J a > - (a 2 + (2a - 6) 2 ) - - (6 2 + (26 - c) 2 ) for all a,b,c £ R, 

y — 2/2 2 

which implies the G-stability of the BDF2 method; see [6] and Lemma The entropies 
i? Q and are formally related by E^[nk, rik-i] = E a [rik\ + 0(r) as r — >■ for k > 2. 

Corollary 2 (Entropy dissipation). Let the assumptions of Theorem U\ hold for a > 1. 
TTien t/ie scheme is entropy dissipative in the sense of 

(13) E%[n k+1 ,n k ] + 2n a r [ (A^jfdz < ^[n fc ,n fc _i], fc > 1. 

In particular, k H- E^[rik, rik-i] is nonincreasing. 

We stress the fact that the implicit Euler scheme ([2]) dissipates all admissible entropies, 
whereas the BDF2 scheme just dissipates one entropy, E^[nk\, where a has been fixed in 
the scheme. 

The proof of Theorem [1] is based on the semi-discrete entropy stability inequality (1T2"]) 
and the Leray-Schauder fixed-point theorem. Instead of (J7J), we employ the algebraic in- 
equalities (fl8l) and ( |T9|) (see Section [2]). We have not been able to obtain similar inequalities 
for BDFfc methods with 3 < k < 6. The reason might be the fact that the only G-stable 
BDF methods are the BDF1 (implicit Euler) and BDF2 discretizations [6]. Moreover, we 
have not been able to prove entropy dissipation for a = 1 since in this case, inequalities 
( TT5|) and ( TT§|) cannot be used. 

If a = 1, we prove that the semi-discrete solution to the BDF2 scheme converges to the 
continuous solution with second-order rate. 

Theorem 3 (Second-order convergence). Let the assumptions of TheoremUlhold, let a = 1, 
and let (vk) be the sequence of solutions to (I51)- ([TD1) constructed in Theorem^ We assume 
that there exist values fik > such that > \x\. > in T d . Furthermore, let n be a strictly 
positive solution to © satisfying y/n G H 3 (0, T; L 2 (T d )) D W 2 >°°(0, T; L 2 (T d )). Then there 
exists a constant C > 0, depending only on the L 2 (0, T; L 2 (T d )) norm of (y/n) tu , the 
L°°(0, T; L 2 (T d )) norm of (y/n) tt , andT, but not on r, such that 



IK - y/n(t k , Oliver*) < Cr 2 , 
where 0<r< 1/8 is the time step and tj, = rk, k > 0. 

It is shown in [21 Theorem 6.2] that the solution n to (j2J) is smooth locally in time if 
the initial datum is positive and an element of W 1,oc (T d ). The proof of Theorem [3] is 
based on local truncation error estimates and the monotonicity of the formal operator 
A(v) = v 1 ~ 2 / a d 2 j(v 2 d 2 j logf) for a = 1 [26]. If a ^ 1, the operator A seems to be not 
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monotone, and our proof does not apply. Possibly, the second-order convergence for q^I 
could be achieved by applying suitable nonlinear semigroup estimates. 

Next, we investigate a fully discrete numerical scheme which dissipates the Fisher infor- 
mation. To this end, we employ the discrete variational derivative method of Furihata and 
Matsuo [13]. The method is based on the variational structure of the DLSS equation, 



The dissipation of the Fisher information F[n] (see 01])) follows from (formally) integrating 
by parts in 



(see (jSJ)). The idea of the method is to derive a discrete formula for the variational derivative 
5F[n]/5n in such a way that the above integration by parts formula and consequently the 
dissipation property hold on a discrete level. We provide such formulas for spatial finite 
difference and temporally higher-order BDF approximations. 

The numerical approximation for equation (TjQ), derived in [9], takes advantage of the 
gradient-flow structure in the sense that the variational structure was discretized instead 
of equation ([1]) itself. The method is based on the minimizing movement (steepest descent) 
scheme and consequently dissipates the discrete Fisher information. In each time step, a 
constrained quadratic optimization problem for the Fisher information needs to be solved 
on a finite-dimensional space. Each subproblem has to be solved iteratively, leading to a 
sequential quadratic programming method. In general, this structure-preserving approach, 
known as "first discretize, then minimize" , has good stability properties and captures well 
other structural features of equations, like those presented in |29j . 

The strategy of the discrete variational derivative method is the standard "first mini- 
mize, then discretize" approach, i.e., the discretization of equation (pQ), as the minimality 
condition in the variational setting, is performed. To some extent this is simpler than the 
above approach, since in each time step only a discrete nonlinear system has to be solved, 
and the main structural property remains preserved. Furthermore, we derive temporally 
higher-order discretizations, whereas the scheme in [9] is of first order only. 

To simplify the notation, we consider the spatially one-dimensional case only. The 
extension to the multidimensional situation is straightforward if we assume rectangular 
grids. Let x , . . . , x^ be equidistant grid points of T with mesh size h > and x = xjy. Let 
U]° be an approximation of n(t k , Xj) and set U k = (Uq, . . . , C7"Jy-_ 1 ) , Un = Uq. Furthermore, 
let Si' 9 be the g-step BDF operator at time for instance, 



(14) 




t > 0. 




(16) 



(15) 




APPROXIMATIONS OF A QUANTUM DIFFUSION EQUATION 



7 



We denote by 8\ the central finite-difference operator at x iy i.e. 8^ U k = (U^ +1 — U£_^)/h. 
Then, following (Tl4|) . we propose the fully discrete scheme 

(17) S^Ut» = 6? SFd tUk ^ l) )) ■ *>*-!, 

where z = 0, . . . , N - 1. The discrete variational derivative 8F d /8(U k+1 , . . . , e 
is defined in such a way that a discrete chain rule holds (see (EHD and (1121) in Section [3] for 
the precise definitions), yielding the dissipation of the discrete Fisher information F d [U k ] 
in the sense of the following theorem. 

Theorem 4 (Dissipation of the Fisher information). Let N G N, U° G WL N be some 
nonnegative initial datum with unit mass, 

Eio u i h = 1 > and letU 1 ,..., U q ~ l G R N be 
starting values with unit mass and i^ff/ 9 " 1 ] < • • • < F d [U°] < oo. Then scheme (fT7|) . with 
the discrete variational derivative 8F^/8{U +1 , . . . , U k ~ q+1 ) defined by ( I42p . consistent 
of order (q, 2) wzt/i respect to the time-space discretization. Furthermore, U k is bounded 
uniformly in k, has unit mass, and the discrete Fisher information is dissipated in the 
sense of 

8l' q F d [U k ] < for all k > q. 

Furthermore, forq = 1 the discrete variational derivative is defined by fl39|) . scheme (fTTj) is 
consistent of order (1, 2) and the discrete Fisher information is nonincreasing, Fd[U k+1 ] < 
F d [U k ] for all k > 1. 

We say that a scheme is consistent of order (q, m) if the truncation error is of the order 
0(r q ) + 0(h m ) for r ^ and h 0. 

The paper is organized as follows. The analysis of the BDF2 time approximation is 
performed in Section [2], and Theorems [1] and [3] are proved. The fully discrete variational 
derivative method is detailed in Section [31 and Theorem |4] is proved. Numerical experi- 
ments in Section H] illustrate the entropy stability, entropy dissipation, and energy (Fisher 
information) dissipation, even in situations not covered by the above theorems. 

2. BDF2 TIME APPROXIMATION 

First, we collect some auxiliary results. The following lemma is needed to show a priori 
bounds for the semi-discrete solutions to the DLSS equation. 

Lemma 5. It holds for all a, b, c G M, 

(18) 2 (-a - 2b + -c) a > -a 2 - 2b 2 + -c 2 + (a - b) 2 - (b - c) 2 , 

(19) 2 - 2b + ~c) a > 1 (a 2 + (2a - 6) 2 ) - ~ (6 2 + (26 - c) 2 ) . 
Proof. We calculate 

2 (-a - 2b + \c ) a = -a 2 - 2b 2 + -c 2 + (a - 6) 2 - (6 - c) 2 + -(a - 26 + c) 2 , 
V 2 2/2 2 2 
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which proves the first assertion. Because of 

2 Qa - 2b + a = ± (a 2 + (2a - bf) - l -{b 2 + (26 - cf) + I( a - 26 + c) 2 , 

the second assertion follows as well. □ 

We also recall the following inequality (see [221 Lemma 2.2] for a proof). 

Lemma 6. Let d > 2 and y/n G H 2 (T d ) n W 1A (T d ) n L°°(T d ) urctfc inf T d n > 0. Then, for 
any (Vd - l) 2 /{d + 2) < a < (y/d + l) 2 /(d + 2), a ^ 1, 

77 1 n / nd^logn^n^dx > K a / (An a / 2 ) 2 dx 
4(a - 1) Jjd Jjd 

and for a = 1, 

- I n(<9 2 (logn)) 2 dx > m / (Av^) 2 dx, 

w;/jere 

P(a) , / n 2 2(d+l) /d- l x 
K a = — — > and p(a) = —a H — —a 



a 2 (p(a) -p(0)) ^ v ' d + 2 V rf + 2 

Proof of TheoremUi Given vq = n^ 2 , the existence of a nonnegative weak solution v\ G 
# 2 (T d ) to (HUD is shown in [22]. Assume that v 2 ,...,v k e H 2 (T d ) are solutions to ([IT]). We 
introduce the variable y by f^+i = e ay l 2 such that n^+i = e y . First, we prove the existence 
of a weak solution y G i/ 2 (T d ) to the regularized equation 

(20) ^ e(W% (l eay/2 - 2 ^ + + l d U eVd iv) + eL (y) = o. 

where e > and 

= A 2 y-div(|Vy| 2 Vy) + y. 

1: Definition of the fixed-point operator. Given z G W 1)i {T d ) and a G [0,1], we 
define on H 2 (T d ) the forms 

a(y, 4>) = \ I e z d 2 3 yd 2 ^dx + e I (AyA0 + | Vz\ 2 Vy ■ V0 + #)dx, 

m = -—( e^ 2 > (\e^ 2 - 2v k + 1^ 0dx. 

Since H 2 (T d ) W 1,4 (T d ) L°°(T d ) with continuous embeddings (remember that d < 3), 
these mappings are well defined and continuous. Furthermore, by the Poincare inequality 
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for periodic functions with constant C P > 0, the bilinear form a is coercive, 

4v\\h(W = e [ (\V 2 y\ 2 + |Vy| 2 + y 2 )dx < e [ {{C 2 P + l)|V 2 y| 2 + y 2 )dx 

JT d J7 d 

= e [ ((C 2 + l)(Ay) 2 + y 2 )dx < e(C 2 P + 1) / ((Ay) 2 + y 2 )dx 

JT d JT d 

<(C 2 P + l)a(y,y). 

By Lax-Milgram's lemma, there exists a unique solution y £ H 2 (T d ) to 

a(y,0) = /(0) forall0£ tf 2 (T d ). 

This defines the fixed-point operator S : W 1A (T d ) x [0,1] -> W^T*), a) = y. It 
holds S'd/, 0) = for all y £ W 1,A (T d ), and 5 is continuous and compact, in view of the 
compact embedding H 2 (T d ) W 1,4 (T d ). In order to apply the Leray-Schauder theorem, 
it remains to show that there exists a uniform bound in W 1,A (T d ) for all fixed points of 
S(;a). 

Step 2: A priori bound. Let y £ H 2 (T d ) be a fixed point of S (•,<?) for some cr £ [0, 1]. 
We employ the test function = y in (1201) . This gives 

= % L e(W% (r Q " /2 _ 2Vk + 1^- 1 ) ydx 

(21) + \f ey(dly) 2 dx + e [ ((Ay) 2 + | Vy| 4 + y 2 )dx. 

To estimate the first integral, we distinguish the domains {y < 0} and {y > 0}: 

— / (\eyi 2 - 2v k + U^) y dx 

ar J Jd \2 2 J 

= — I (3e y y - Ae {1 ~ a/2)y v k y + e {l - a/2)y v k ^y)dx 
aT J{y<o} 

+ — I (3e y y - Ae {1 - a/2)y v k y + e {1 - a/2)y v k ^y)dx. 

& T Jiv>0\ 



The first integral on the right-hand side is estimated by using the Young inequalities 

_ 4e (l-a/2)y Vky > _ 2e (2-a)y y 2 _ 2y 2 and e (l-a/2)^ fc _ i?/ > -^-^Vy^ _ 



— [ (3e y y - 4e {1 ~ a/2)y v k y + e {1 ~ a/2)y v k ^y)dx 

aT Jiv<0\ 

J <q fey - \e^ y y 2 - 2v\ - \v 2 _^ dx 



'{y<o} 

a 

> — 

ar 



e y (y - 1) + 1 + (1 + 2y) e y - |e (2 " Q) y -\-2v\- \v\_ x ) dx. 



ar J {y<0} \ 
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Since y h->- (1 + 2y)e y — ^ 2 ~ a ^ y y 2 — 1, y < 0, is bounded from below (remember that 
a < 2), we find that 

— / (3e y y - 4e {l - a/2)y v k y + e {1 - a/2)y v k . iy )dx 

[ {e y (y - 1) + l)dx - — Cl - — / Uv\ + \v\\ dx, 

J{y<0} « r « r J{y<0\ V 2 / 



'{s/<o} 

a 

> — 

ar 



where ci > depends only on the lower bound of y t— >• (1 + 2y)e y — |e^ 2 "^y 2 — 1, y < 0, 
and meas (T d ). For the remaining integral over {y > 0}, we employ the Young inequalities 
_ 4e (i-a/2)2/^ y > _ 2e (2-*)y -y^-vl and e^-^v^y > -\e^-^ y - ±?/ 4 - \v\_ x : 



— I (3e y y - Ae {1 - a/2)y v k y + e {1 - a/2)y v k . iy )dx 
ar -'{y>o} 



> — / f 3e»j/ - -e( 2 "^ - -y 4 - ^ 4 - -vt ,) dx 
-arrj {v > 0} \ 2 4 y fe 4 fc "V 

/ (e y (y - 1) + 1 + ((1 + 2y)e* - h^ y - V - l) 



a 

ar 



^ - T^-i) da; - 



The mapping y i— >• (1 + 2y)e y — le 1 - 2 a ^ — |?/ 4 — 1, y > 0, is bounded from below which 
implies the existence of a constant c 2 > such that 



a 

ar 



(3e y y - 4e {1 - a/2)y v k y + e {1 ~ a/2)y v k _ iy )dx 
{y>o} 



a 

> — 

ar 



/ (e y (y - 1) + l)dx - — c 2 - — / L 4 + dx. 

J{ y >o} ' ar ar J {y > 0} \ 4 J 

Summarizing the estimates for both integrals over {y > 0} and {y > 0}, it follows that 

(22) — / e^^ y (%^ 2 - 2v k + ydx > — [ {e y (y - 1) + l)dx 

ar J T d \2 2 J ar J T d 

- — [ \2v 2 k + v% + -vl_ x + -vl_A dx - — (ci + c 2 ). 
ar 7 T d V 2 4 y ar v y 

For the second integral in (|2~TT) . we use Lemma [6j 

1 f e y (d 2 3 y) 2 dx > 2m [ (Ae y/2 ) 2 dx, 

2 J^d Jjd 
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where k± > depends only on the space dimension d. With this estimate and f l22|) . equation 
($IT§ implies that 

— / (e y (y-l) + l)dx + 2K 1 f (Ae y/2 ) 2 dx + e [ ((Ay) 2 + | Vy| 4 + y 2 )dx 

(XT J T d Jjd Jjd 

< — I [2v 2 k + vt + -vl_i + T Vfc_i J dx + — (ci + c 2 ). 

By the definition of the entropy, this inequality can be written as 

El [n} + ^^ [ (A e y/ 2 ) 2 dx + — [ {(Ay) 2 + \Vy\ 4 + y 2 )dx 
a J T d a J T d 

(23) < [2vl + vt + ^v 2 k _ x + l -v A k _^j dx + c 1 + c 2 . 

The right-hand side gives a uniform (with respect to a) bound since i>fc_i, v k G W 1,A (T d ). 
Hence, by the Poincare inequality we obtain the i/ 2 -bound 



\y\\ 2 H 2 iTd) < C [ ((Ay) 2 + y 2 )dx<C, 



where the constant C > depends on a, e, r, v^, and but not on a. The continuous 
embedding H 2 (T d ) c — >■ iy 1,4 (T d ) then implies the desired uniform bound, w rl - 4 (Tr rf ) < C- 
Leray-Schauder's fixed-point theorem provides the existence of a fixed point y e of S^y, 1) = 
y, i.e. of a solution to fTZOj) . 

Sfep 3: Limit e — > 0. Let ?/ £ be a solution to ( |20|) . constructed in the previous steps. Set 
w e := e aye l 2 and n e := e Ve . Then t> £ solves 
(24) 

l/J*- 1 Q« e - 2 Wfc + + flj (^e 27 ^ 1 ^.^ - ^(^^A*)) + eL(y £ ) = 0. 

The goal is to pass to the limit e — > in this equation. 

Let a > 1. We employ the test function e (Q - 1)y V(a — 1) G H 2 (T d ) in (|2QJ) and find that 

0= , 2 n / f - 2^ fe + ^ £ dx + — / e^^e^^d* 

a(a - 1) i T d V 2 2 / 2(a-l)y T d J 7 

+ ^-(L( 2/£ ),e^- 1 ^)^ 2iH2 . 
a — 1 

Inequality f TTHT) shows that 

2 /" /3 1 \ 1 /" /3 2 2 1 2 \ 

— — / -u £ - 2w fc + f £ dx > — — / l-v E - 2v k + -t» fe _ 1 dx 

a(a-l)J T d\2 2 J a(a - 1) J T d \2 £ 2 \) 

+ p 1 / ((u e - ^fc) 2 - K - ^fc-i) 2 )dx. 
a(a - 1) J T d 

The integral involving the second derivatives is again estimated by using Lemma [6j 

T [ e^dly £ d 2 Ae^- l ^)dx>2K a r [ (Ae^' 2 ) 2 dx = 2n a r [ (Av £ ) 2 dx. 
2(a-l) J T d J T d J T d 
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Now let us consider the e-term and show that (L(y £ ), eS a ~ x ' Ve )H- 2 ,H i is bounded from 
below uniformly in e. By construction, v £ and n £ are strictly positive since y £ G H 2 (T d ) ^-t 
L°°(T d ). Therefore, we can write (cf. [221 Section 4.1]) 



(L(y £ ), e^^H-tjp = 4(a - 1) / (e^ 1 * ( 

j T d \ \ 



-(2 -a) 



2\ 2 



dx 



+ 4(a 2 - 1)(3 - a) [ e^" 1 ^ 4 dx + / y e e (Q - 1)!/£ dx > -C, 

Jjd e^l 1 Jjd 

where C > depends only on a. We have used the fact that le'"" 1 ' 1 > — 1/ ((« — l)e) for 
all ieK. 

Summarizing the above inequalities, we obtain 

-7- 1 n / - 2 ^ 2 + ^Li J d * + zt^t / ( ( ^ _ Wfc)2 " _ ^-i) 2 ) dx 

(25) + 2rn a \ (Av £ ) 2 dx < Ce. 



Inequality (125j) provides the estimate for (v £ ) in H 2 (T d ) uniformly in e. Therefore, there 
exists a limit function v G H 2 (Y d ) such that, up to a subsequence, as e — > 0, 

w e ^t; weakly in H 2 (T d ), 

v £ ^v strongly in W lA (Y d ) and L°°{T d ). 

Consequently, since 2/a — 1 > 0, 

(26) v 2 J a - 1 d 2 ij v £ -± v 2ta ~ x d 2 jV weakly in L 2 (T d ), i,j = l,...,d. 

According to the Lions- Villani lemma on the regularity of the square root of Sobolev 
functions (see the version in [31 Lemma 26]), there exists C > independent of e such that 

II VVeWw 1 ' 4 ^) — C|| u e||_ ff 2 (T<*) < C. 

Since 1/2 < 1/a < 1, Proposition A.l in [23] shows that the strong convergence v £ — > t> in 
//' 1 (T d ) and the boundedness of (y'tQ in W 1,4 (T d ) imply that 

v l/a _^ v l/a strongly in W l,2a(jdy 

Hence, we have 

(27) diivl'^d^vl 101 ) -> ^(w 1/Q )^(w 1/a ) strongly in L a (T d ), i, j = 1, . . . , d. 
Estimate (123]) and Ei[n] > provide the uniform bound 

y/£\\ye\\H*(T*) + V^IIVt/eH^^d) < C, 

which shows that 

(28) eL(y £ )^0 weakly in £r 2 (T d ). 
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Using (j) G W 2,OQ (Q) as a test function in the weak formulation of (|24p . the convergence 
results fEo^ -f l^S]) allow us to pass to the limit e — > in the resulting equation, which yields 
fTTT]) for Vk+i ■= v . In fact, it is sufficient to use test functions <p £ W 2,a ^ a ~ l " > {T d ). 

If a = 1, the convergence result follows similarly as above based on the uniform bound 
||e y£ / 2 ||#2 < C, which is obtained from a priori estimate ( |23|) . using the elementary inequal- 
ity s < s(logs — 1) + e for all s > 0, which gives a uniform L 2 -bound for e Ve . In that case, 
the test functions <fi G H 2 (T d ) can be used in 

Step 4 : Entropy stability. Let a > 1. Using the test function v\ 2 ^ a /(a — 1) in (ITUl) . it 
follows that 



Ta(a — 1) Jjd 2{a — 1) J T d J J 

By Lemma [6], we infer that 

(29) 1 / (^ + (^- V) 2 )dx + 2™ Q / (A Vl ) 2 dx < [ v 2 dx. 

This gives an i7 2 -bound for t>i. 

Next, let k > 1 and let y s be a weak solution to (1201) . Set t> £ = e ay ^ 2 . The convergence 
results of Step 3 allow us to pass to the limit e — > in (1231) . Using the weakly lower 
semi-continuity of u H- ||Aw|| 2 2 ( T d) on H 2 (T d ), it follows that 



a (a — 1) Jjd \2 2 



+ , 1 7T / ((ufc+i - - K - ffc-i) 2 )dx 
a(a - 1) J T d 



(30) + 2k q t / (Av k+1 ydx < 0, 

where, as before, v k +i = liTa E ^.oV E . Summing (129]) and (130|) over k = l,...,m — 1, some 
terms cancel and we end up with 

m— 1 



2a(a 



- f v 2 m dx + lK a r V [ (&v k+1 ) 2 dx < / f + + ^ 2 )dx. 

1) J T d J T d 2a{a - 1) i T d 



Set a m = ||"ym|| 2 2( n ) for m > 0. By (1291) . ai < a . Then, the above estimate shows that 
&m < l^m-i + §Qo- A simple induction argument gives a m < ao for all m > 1. Therefore, 

~ 77 / ^i + JvV / (Av k ) 2 dx < f v 2 dx. 

This implies the entropy stability estimate (fl2|) . □ 
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The proof of Corollary [2] is a consequence of the above proof. Indeed, employing in- 
equality ffTUj) instead of fTTH]) . we can replace f[3TIj) by 

7r~T~ — 7T / ( v l+i + fak+i - v k ) 2 )dx + 2n a r [ (Av k+ i) 2 dx 
2a(a-l)J T d x ^ J T d 

which equals (fT3|) . 

Next, we prove that, if a = 1, the solutions v k are smooth as long as they are strictly 
positive. 

Lemma 7. Let a = 1 and /et (?;&) be the sequence of weak solutions constructed in Theo- 
rem[l\ satisfying Vk > fi k > in T d for k > 1 and some fi k > 0. T/ien v*. G C°°(T d ). 

Proof. We recall that the weak form ( Hip for a = 1 reads as 

Vk+i Q^fc+i _ 2t; fc + Tj^fc-iJ da; + ^ y (^fc+i^^fc+i - diVk+tdjVk+^dPjCfxlx = 
for G if 2 (T d ). Since v k is assumed to be strictly positive, we can write 



1 

2 

,2 



Ufc+i^fc+i - diV k+1 djV k+1 = -ra fc+ i^ 7 logn fc+ i, 



(31) -v fc+1 - 2u fc + -Ufc.! + -dUrik+id*, \ogn k+1 ) = in #~ 2 (T d ). 



where n fc+1 = and consequently, 

' 3 1 ^ 

-v k+1 - 2v k + -Vk-! 1 + -c vV , K ^ v 

With the identity 

dfjUk+idjUk+t (d j n k+1 ) 2 d i n k+1 



5y(n fe+1 ^ logn fc+ i) = A ra fe+1 - d t 2- 

v n i;+i 

it follows that n k+ \ solves 

(32) A n fc+1 = 2^ J - - -v k+1 -v k+1 - 2v k + -v k _x 

\ n k+1 n 2 k+1 ) t \2 2 

in the sense of H~ 2 (Y d ). The second term on the right-hand side is an element of H 2 (Y d ). 
The continuity of the Sobolev embedding H 2 (T d ) <— > W l,& (Y d ) (for d < 3) implies that 
(d jnk+1 ) 2 d ink+1 /n k+1 G L 2 (T d ) and ^.n fc+1 ^n fc+1 /n fc+1 G L 3 / 2 (T d ) H' 1 / 2 ^) for all 
i,j = l,...,d. This proves that 

A 2 n k+l G ir 3/2 (T d ). 

The regularity theory for elliptic operator on T d (e.g., using Fourier transforms on the 
torus) yields n k+ i G H 5 ^ 2 (T d ) which improves the previous regularity n k+ i G H 2 (T d ). 
Taking into account the improved regularity and the embedding H 5 ^ 2 (T d ) <-> W 2 ' 3 (T d ), we 
infer that the right-hand side of ( 1321) lies in i? _1 (T d ), i.e. 

A 2 n k+1 eH~\T d ), 
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which implies that Tik+i G H 3 (T d ). By bootstrapping, we conclude that n k+ i e H m (T d ) 
for all m G N. □ 

Now, we are in the position to prove Theorem [3J 

Proof of Theorem^ Let (v k ) be a sequence of weak solutions to ( ITTj) . Since we have 
assumed that is strictly positive, Lemma [7] shows that v k is smooth. As a consequence, 
Wfc solves (see (151) ) 

3 11 

-v fc+1 - 2w fc + -w fc _i + <9j log u fc+ i) =0 in T d . 

Let n = v 2 be a solution to (j5J) with the regularity indicated in the theorem. By Taylor 
expansion, 

1 /3 1 \ f 

v t (t k+ i) = - f g «(*fc+i) " 2 «(*fc) + Tf(^-i) ) + 7' fc ^ 



where 



fk = - vttt(s)(t k - s) ds + - / Uttt(s)(*fc-i - s) ds 
can be interpreted as the local truncation error. We estimate f k as follows: 

m— 1 

(33) 2^ ll/fcllL2(Td) — C^II' i; ^lli2( 0iT;L 2( T[ i))r 5 , 



where Cr > does not depend on r or m. Similarly, we have 

Vt(h) = -(v(t 1 )—v(t )) + —, where / = / %(s)sds, 
r r 



and 



T 2 



(34) ||/o|U2( T d) < y ||f tt (s)|| i 2 (T d ) sds < y\\ v i*\\l°°(o,T;V(j«))- 

Replacing the time derivative v t in ([2]), written as t>j + t> ~ 1 9 J 2 J (f 2 <9 2 - logv ) = 0, by the above 
expansions, it follows that 

(35) v(t 1 )-v(t ) + ^d 2 ] (v(t 1 ) 2 d 2 J \ogv(t 1 )) = -f , 

(36) - 2v(t k ) + \v{t k ^) + — — <9 2 {v{t k+l ) 2 dl \ogv{t k+l )) = -f k , 

for k > 1. Taking the difference of ( fTOl) . multiplied by v± , and ( l35l) . and the difference of 
dHJ), multiplied by and (1361) . we obtain the error equations for e k := v k — v(t k ): 

ei - e + r(A(ui) - = /o, 

3 1 

-e k+1 - 2e k + -e fc _i + r(A(w fc+1 ) - A(u(£ fe+1 ))) = / fc , A; > 1, 
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where we have introduced the operator 

A : D(A) H~ 2 (T d ), A(v) = -dWd^ logv), 

v 

with domain D(A) = {v G H 2 (T d ) : v > in T d }. 

We multiply the error equations by e\ and e k +i, respectively, integrate over T d , and sum 
over k = 0, . . . , m — 1: 

/ (ei - e )eidx + V" / (\e k+1 - 2e k + -< 
(37) + r / W u fc+0 ~~ ~ ^(4+i))dx = ^ / fc e fc+ ida;. 



-e fc _i ) e fc+ idx 



m— 1 



fc=0 



fc=0 



Using e = and inequality ffTSl) . the first two integrands can be estimated by 

m-i ,^ 1 v 

(ei - e )ei+ f -e fc+ i - 2e fc + -e fc _i J e fc+ i 
k=i ^ ' 

m ~ 1 /3 1 1 1 

- e i + ( i e fe+i _ e '- + i e '~i + 2 (6fc+1 ~ 6k)2 " 2 ' 
fc=i ^ 



e k - e fc _i) 2 ) 
1 



= e? + -e& - |je 2 - JeJ^ + ^-e 2 , + ^(e m - e m _i) 2 - ~(ei - e ) 2 

~~ ^ m ^ c m-l 4 1 2 °m-l) 
> L 2 _ l e 2 _ 1p 2 

- 4 m 4 4 !• 

For the third integral in (l3Tj) . we employ the monotonicity of the operator A. In fact, it is 
proved in [26j Lemma 3.5] that for positive functions Wi, w 2 G i/ 4 (T d ), 



(A(ioi) - A(w 2 ))(w 1 - w 2 )dx 



div wfV 



The right-hand side of (137]) is estimated by Young's inequality: 



w 1 - w 2 

Wx 



dx > 0. 



/ /o^ldx < 2||/o||| 2 fT^ + o ll e lllL 2 (T d )' 
J T d 8 

^ fk^k+idx < —WfkW 2 ^^ + -||efc + i|| 2 2 ( T d), 



k > 1. 
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Summarizing the above estimates and taking into account fl33l and (|34p . we find that 
3 11 1 

ll li'"* ll 1 1 2 ll li'"* "i 1 1 "* 1 1 '"* ll 1 1 2 

^||e m || L 2 (T d) < ^||e m --i|| L 2 {T d ) + ^||ei|| L 2( X d) + 2||/ || L 2( T d) + - ||ei|| L 2 (T d) 

^ m— 1 m— 1 

+ 7^1 H/ fe Hi 2 (T d ) + 2 ll e fc+ 1 lli a (T <1 ) 

fe=l k=l 

1 3 m 

In i|2 ^ii i|2 ^~ \ ^ II II 2 

^ -|H e m-l lli 2 (T d ) + g ll e lllL 2 (T d ) + ^ T + 77 H efc llL 2 (T d )' 

fc=2 

where C > depends on the L 2 (0,T; L 2 (Y d )) norm of ^ and the L°°(0,T; L 2 (Y d )) norm 
of Vtt but not on r. Taking the maximum over m = 1, . . . , M, we infer that 

3 1 r A 

Il2 ^ ^ i|2 /^y 4 X II Il2 

7 maX », H em nL 2 (T d ) ^ o max . f ll e m-l|lL 2 (T d ) + ^ T + p / , II e fc II L 2 (T d )" 

4 m=l,...,M v ' o m=l,...,M v ' Z ^ 

fc=2 

The first term on the right-hand side is controlled by the left-hand side, leading to 

M 

|| e M 1 1 Z.a/Tpd-i < max || e m|li2("]rd-) < 8CV +4r\ l^fclli^lW 
v ' m=l,...,M v ' ^— ' v 1 

k=2 

We separate the last summand in the sum, 

M-l 

(1 - 4r)||e M ||^ (Td) < 8CV 4 + 4r ^ |N|^ 2(T(i) , 

fc=2 

and apply the inequality 1 + x < e x for all x > and the discrete Gronwall lemma (see, 
e.g., [301 Theorem 4]): 



2 8C7r 4 / 4r \ M ' 2 8CV 4 / 4^ x 4 



The result follows for all < r < 1/8 with the constant AyC exp(4T), where T > is the 
terminal time. □ 



3. Fully discrete variational derivative method 

In this section, we explore the variational structure of the DLSS equation on a discrete 
level, using the discrete variational derivative method of [13] . In order to explain the idea, 
we consider first the implicit Euler discretization. 

Let Xi — ih, i — 0, . . . , N — 1, be an equidistant grid on the one-dimensional torus T = 
[0, 1), let t k = kr with r > 0, and let approximate n(t k , xi). Set U k = (17* , . . . , f/jv-i) e 
WL N and Ug = Ug mo d n for all £ G Z. We introduce the following difference operators for 
U = (Ui) G R N : 
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forward difference: 8fU = h~ x (U i+ i — C/j), 

backward difference: 5~U = — U^i), 

central difference: 5^U = (2h)- 1 (U i+1 - U^), 

second-order central difference: 5f U = S^5^U = 5~5fU. 

The first step is to define the discrete Fisher information. We choose a symmetric form 
for the derivative, v 2 x {xi) « \{{5fV) 2 + (S^V) 2 ), where V = (V<) = (V^i) e The 
Fisher information F[w 2 ] = J* T t^dx is approximated by using the first-order quadrature 

rule J T w(x)dx ~ Yl!i=o w i x i)^ 1 - Actually, this rule is of second order 0{h 2 ) here, since 
due to the periodic boundary conditions, it coincides with the trapezoidal rule, (w(x ) + 
w(x N ))h/2 + J2j 

=1 w(xi)h. Therefore, the discrete Fisher information reads as 

N-l 

F AU] = 2 E ^ 5 tvf + (s;v) 2 )h, u = m e r n . 

i=0 

The second step is the definition of the discrete variational derivative. Applying the dis- 
crete variation procedure and using summation by parts (see [131 Prop. 3.2]), we calculate 

N-l 

' h 



F d p k+i ] - F d [u k ] = \J2 (( 5 t yk+1 ) 2 - ( 5 t yk ) 2 + (^y k+1 ) 2 - (^v h ) 2 ) 

i=0 
N-l 

= -J2 W( yk+1 + V k )5f{V k+1 - V k ) 



2 

i=0 
N-l 

2 

i=0 

+ 5'(V k+1 + V k )5-(V h+1 - V k )\h 



N-l 

= ~ Yl 6 ?\ vk+1 + V k ){V k+1 - V k )h 

i=0 

(38) =-Ete?r-* *>o. 

This motivates the definition of the discrete variational derivative 

(39) 5Fd S?\v k+1 + V k ) 2 = Q N _, 

since this implies the discrete chain rule 

F d [U^} - F d [U k ] = £ 6{u L( Uk)i (U k+1 - U k )h. 



i=0 



Observe that (139]) is a Crank-Nicolson type approximation of the variational derivative 
5F[n}/5n = — (\/n) xx / \/n = —v xx /v, where n = v 2 . The implicit Euler discrete variational 
derivative (DVD) method for the DLSS equation is then given by the nonlinear system 
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1!) 



with unknowns U k+1 = (V k+1 ) 2 : 

(40) - U{ ) = (u M S? ( HU S J* ut) ) ) • i = » ~ 1. * > 0. 

The initial condition no is approximated by its projection on the discrete grid, defining the 
starting vector U° G M. N . Multiplying the above scheme by 5Fd/5(U k+1 ,U k )i, summing 
over i = 0, . . . , N — 1, and employing the discrete chain rule (|38p . we infer the discrete 
dissipation property 



In fact, this proves the monotonicity of the discrete Fisher information for 5=1. 

Remark 8. Observe that we could have taken a different approximation for the discrete 
Fisher information, e.g. Fd[U] = ^i^i^Vyh. This would lead to a different variational 
derivative SF d /5(U k+1 ,U k ) and eventually to a another scheme (HOjl . with F d replaced by 
Ffa which dissipates F^ instead. Besides the symmetry, which brings the second-order 
consistency in space, the above choice of the discrete Fisher information is motivated by 
the fact that 8f8^ = df \ used in the discrete variation procedure. □ 

In the following, we consider temporally higher-order discretizations. There are several 
ways to generalize the above DVD method. In order to stay in the spirit of Section |2} we de- 
rive higher-order DVD methods, which are based on backward differentiation formulas. The 
function /(£, rj) = (£ 2 + rf)/2 represents both the Fisher information F[n] = f T f(v x , v x )dx 

and the discrete Fisher information Fd[U] = Yl!i=o fi^t^^i^)^ 1 - The definition of / is 
motivated by the following formal representation of the variational derivative, 



This formula gives an idea how to approximate the variational derivative in general. We 
denote by Sl' q the q-ih step BDF operator at time t^. For instance, the formulas for q = 1 
and q = 2 are given in ffl5l) and ffl6]) . respectively. The discrete variational derivative of 
order q is defined componentwise by 

( 42 ) = {S ^ f) + St{ ^ f)) ' k ~ q ~ 1j 

where the discrete operators d^f and &zf are given by 

\^=8+V k + 
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and r corr is a correction term, which has to be determined in such a way that the discrete 
chain rule 

holds. The role of the correction term is not only to satisfy the discrete chain rule but also 
to increase the temporal accuracy of the discrete variational derivative. Straightforward 
computations with the above expressions using summation by parts formulas and periodic 
boundary conditions yield 

5F d _ Sfv^ jfojfV^ t> 

5{U k +\ U k -^ l )i ~ V k+l Tcorr V k+1 ' q 



We note that for q = 1, this definition generally does not coincide with the discrete varia- 
tional derivative ( 1391) . The temporally BDFg discrete variational derivative (BDFg DVD) 
method is then defined by the following nonlinear system in the unknowns U k+1 = (V k+1 ) 2 : 

(45) S&VT 1 - 4P> (p»»f, W ( f (gtH[ m ut - q+1) )) ■ i = 0,...,W-l, *>,-!. 

In particular, for q — 1, we obtain two methods: the BDF1 DVD scheme (I45p and the 
DVD scheme (HQ]). 

Proof of Theorem^ Let n = v 2 be a smooth positive solution to ([I]) with d = 1. According 
to [2], such a solution exists at least in a small time interval if the initial datum is smooth 
and positive. Furthermore, let q G N, q > 2 (and typically q < 6), be the order of the 
backward differentiation formula. 

First, we consider the discrete variational derivative ( I39p . A Taylor expansion around 
(t k +i,Xi) yields 



5F d 



6l 2) (v(t k+ i,Xi) + v(t k ,Xj)) _ 
6(n{t k+1 ),n(t k )) x= Xi v(t k+1 ,Xi) +v(t k ,Xi) v 

= —[n\(t k+1 ,x, 

where % = 0, . . . , N — 1, fc > 0. Similarly, 



{t k+1 ,Xi)+0(r) + 0(h 2 ) 



SF 

■^[n](t k+1 ,Xi) + 0(T) + 0(h 2 ) 



6 Fa 



5(n(t k+ i),n(t k )) 



Thus, the local truncation error of the right-hand side in (139|) is of order O(r) + 0(h 2 ). 
Since the left-hand side is of order O(r) in time and exact at spatial grid points x», the 
local truncation error of scheme (139|) is of order 0(t) + 0(h 2 ). The monotonicity of the 
discrete Fisher information is shown in (l4Tj) . 
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The mass conservation is an obvious consequence of the scheme. To prove the uniform 
boundedness, we observe that, by the discrete if 1 -seminorm, 



JV-l 



^2{5tV k fh < F d [U°] < oo for aU k > 1. 



Then, according to the discrete Poincare-Wirtinger inequality, for i = 0, . . . ,N — 1, k >1, 
pT3l Lemma 3.3], 

N-l 

\V? - M k \ 2 < J2( 5 tV k ) 2 h < F d [U°] 

with M k = YliLo 1 ^i^ 1 - Jensen's inequality for the quadratic function and the mass con- 
servation property of the method give M k < 1 for all A; > 0. Finally, by the triangle 
inequality, |Vf | < F d [U°] 1/2 + 1 and thus, \U? | < 2F d [U°) + 2. 



Next, we consider scheme ( Hoi) with the discrete variational derivative (I43j) . By construc- 
tion, the left-hand side of ( I43p is of order q in time and exact at the spatial grid points 
Xi. Thus, it remains to prove that the right-hand side is of order (g, 2) with respect to 
time-space discretization. 

Taylor expansions show, with a slight abuse of notation, that 

(46) 5fv(t k+1 , x^ = v x (t k+1 , x^ ± ^v xx (t k+1 , Xi) + 0(/i 2 ), 
47 — = (t fc +i,Xi ) + 0(h , 

V{t k+1 ,Xi) V 
f A o\ ^k+l^i ^(tk+liXi) Tl txx . \ i Q\ i m<l,1\ 

(48) r = (tk+uXi) + 0(T q ) + 0{h ^ 

v{t k+ i,Xi) V 



(49) SfSi^nfa+iiXi) = n tx (t k+1 ,Xi) ± -n txx {t k+1 , x{) + 0(r q ) + 0(h 2 ), 



(50) sf rJkp^ii 

y V(t k+1 ,Xi) 

We prove that r corr is of order (q, 2). Let r n and r d denote the numerator and denominator 
of r corr , respectively, replacing Vt +1 by v(t k +i,Xi) and f/j by n(t k +i,Xi). Taking into 
account the periodic boundary conditions, we find that 



^(<5+c^ 1 n(t fc+1 ,x i ))5, 

i=0 

7V-1 / x l,q 



J2( 6 i KtMtk+U Xi))5. 



v(t k+1 ,Xi) 

^n^+i, Xi 



i=0 



v(t fc+ l,£j) 
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Therefore, we can split into two parts: 

, N-l 



E 

j=0 



v(t k+ i,Xi) 



6 k '^ 1 n(t k+1 , Xi 

v(t k+1 ,Xi) 



In view of (UnD-dSHD, it follows that 



h. 



N-l 



Td 



2j2(nt x v tx ){t k+l ,x i )h + 0{r q ) + 0{h 2 ). 



i=0 



The numerator r n is treated in a similar way. Using fj46l) . the first term in r n can be written 

as 

N-l 



2dt 



AT-l 



i=0 



O(t<0 



k+l 



^(v x v xt )(t k+1 ,Xi)h + 0(r 9 ) + 0(h 2 ). 



i=0 



For the second term in r n , we observe that, because of the periodic boundary conditions, 



Af-l 



^25?v{t k+l ,Xi)5t 



i=0 



v(t k+ i,Xi) 



N-l 



h=^2$i v(t k+1 ,x l )5- 



and hence, employing (146]) and (1501) . 



JV-1 



i=0 



v(t k+1 ,Xi) 



i=0 



N-l 



i=0 



few(tfe+i,a:i) 

f(t fc+ l,Xi) 



5fc+in(tfc+i,Xi) 



+ ^ u(t fc+ i,Xi)5 t : 



U(tjfe+1, Xj) 

2 ^(tttKiw.i^ + O(r') + 0{h 2 ). 



N-l 



i=0 



Summarizing these identities yields r corr = 0{r q ) + 0{h 2 ). Finally, (T47|) -( l48|) imply that 



5F d 



5(n(t k+1 ), . . . ,n(t fc+ i_g)) 



7**i 



SiLdPnitk+^Xi) 



v(t k+1 ,x 
5F[n 
Sn 



v(t k+ x,Xi 

(t k+ i,Xi) + 0(7^ + 0(h 2 ). 



This shows that the discrete variational derivative (1431) is of order q in time, finishing the 
proof. □ 
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4. Numerical examples 

In this section, we present some numerical examples which illustrate the decay properties 
of the entropy functionals and Fisher information as well as the convergence properties of 
the schemes presented in the previous sections. 

4.1. BDF2 finite-difference scheme. The DLSS equation (j2J) is approximated by the 
BDF2 method in time and central finite differences in space. The scheme is given by the 
following nonlinear system with unknowns V k = (U^) a ' 2 : For i = 0, . . . , N — 1 and k — 1, 

{Vt?' a - l (V? - V?) + r8f {{VlfHf log V/) = 

and for i = 0, . . . ,7V- 1, k > 2, 

(lf+ i )2 /.-i (Z_ v h + i _ 2vf + 1^ + T5 f) ((y^fl^ \ogV^) = 0. 

The initial datum (V^°) is given by (no(xi) a ^ 2 ) . For k — 1, the scheme corresponds to the 
implicit Euler discretization, needed to initialize the BDF2 scheme for k > 2. The above 
nonlinear system, with periodic boundary conditions, is solved using the Newton method. 

We choose the initial datum n (x) = 0.001 + cos 16 (7rx), x G [0, 1]. The spatial mesh size 
is h — 0.005 (N = 200) and the time step r = 10~ 6 . The (continuous) entropies E a [n] 
are dissipated for 1 < a < 3/2. Figured] (a) illustrates the stability and, in fact, decay of 
the discrete entropies E a ^, defined below, for various values of a. Although Theorem [1] 
does not provide a stability estimate for a = 1, the numerical results indicate that the 
discrete entropy E\^\JJ] = ^^^(U^logUi — 1) + l)h is decreasing. Figure [1] (b) shows 
that the decay of the discrete relative entropy is exponential, and even the discrete Fisher 
information converges exponentially fast to zero. Here, the discrete relative entropy is 
defined by 

JV-l N-\ 

E r a %U k ] = E aA [U k ] - E a>d [U], where E aA [U k ] = ^(t/f) Q /i, U = ^ U k h. 

i=0 i=0 

According to Theorem [31 the semi-discrete BDF2 scheme converges in second order if 
a — 1. This may be not the case for the fully discrete scheme, since the discretization may 
destroy the monotonicity structure of the spatial operator. However, Figure [2] shows that 
the numerical convergence rate is close to 2, even for q ^ 1. The numerical convergence 
rates cr have been obtained by the linear regression method. The convergence of the 
method is measured in the discrete £ 2 -norm 

(N-l \ X /2 

and the numerical solutions are compared at time t = 5 • 10 -5 . Here, the "exact" solution 
V^i is computed by the above scheme using the very small time step r = 10~ 10 . 
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Figure 1. (a) Entropy stability (decay) for the BDF2 finite-difference 
scheme, (b) Exponential decay of the discrete relative entropy and the dis- 
crete Fisher information for the BDF2 finite-difference scheme. 
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Figure 2. Temporal convergence of the BDF2 finite-difference scheme for 
various values of a; the convergence rate is denoted by cr. 



4.2. Discrete variational derivative method. We present some numerical results ob- 
tained from the DVD and BDFg DVD schemes derived in Section [31 The initial datum 
and the numerical parameters are chosen as in previous subsection. In order to solve the 
discrete nonlinear systems, we employed here the NAG toolbox routine c05nb, which is 
based on a modification of the Powell hybrid method. It turned out that this routine is at 
least three times faster than the standard MATLAB routine f solve. 

In Figure El the temporal evolution of the discrete relative entropies E* [U k ] and the 
discrete Fisher information F^\U ] are depicted for (a) the implicit Euler scheme (140!) and 
(b) the BDF2 scheme (I43j) . We observe that the decay is in all cases exponential. This 
holds also true for the BDF3 scheme (results not shown). 
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Next, we test numerically the convergence in time of the DVD scheme. Figure H] il- 
lustrates the £ 2 -errors of the methods. We have chosen the mesh size h = 0.01, and we 
compared the numerical solutions at time t m = 5-10 -5 . The "exact" solutions are computed 
by the respective method taking the time step r = 10~ 9 . The numerical convergence rates, 
computed by the linear regression method, are given in Table [TJ We note that the BDF3 
DVD scheme gives only slightly better results than the BDF2 DVD scheme. The reason 
is that the first step is initialized by the first-order scheme (140 p . and this initialization 
error cannot be compensated by the higher-order accuracy of the local approximation. In 
order to obtain a third-order scheme, we need to initialize the scheme with a second-order 
discretization. 

10" 3 

10" 4 
=S lO" 5 
= 1CT 6 

lO" 7 

10" 8 

10 -8 1Q -7 1Q -6 

r 

Figure 4. Temporal convergence of the DVD, BDF2 DVD, and BDF3 DVD 
schemes. 
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Scheme 


Convergence rate 


DVD 


1.020 


BDF2 DVD 


1.824 


BDF3 DVD 


1.977 



Table 1 . Numerical temporal convergence rates for the discrete variational 
derivative methods. 
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